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We review some recent developments in numerical algorithms to solve the time-dependent Maxwell 
equations for systems with spatially varying permittivity and permeability. We show that the Suzuki 
product-formula approach can be used to construct a family of unconditionally stable algorithms, 
the conventional Yee algorithm, and two new variants of the Yee algorithm that do not require the 
use of the staggered-in-time grid. We also consider a one-step algorithm, based on the Chebyshev 
polynomial expansion, and compare the computational efficiency of the one-step, the Yee-type, the 
alternating-direction-implicit, and the unconditionally stable algorithms. For applications where 
the long-time behavior is of main interest, we find that the one-step algorithm may be orders of 
magnitude more efficient than present multiple time-step, finite-difference time-domain algorithms. 
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I. INTRODUCTION 



The Maxwell equations describe the evolution of electromagnetic (EM) fields in space and time |jj . They apply to 
a wide range of different physical situations and play an important role in a large number of engineering applications. 
In many cases, numerical methods are required to solve Maxwell's equations E|. A well-known class of algorithms 
is based on a method proposed by Yee Q). This finite-difference time-domain (FDTD) approach owes its popularity 
mainly due to its flexibility and speed while at the same time it is easy to implement [|J |y. 

A limitation of Yee-based FDTD techniques is that their stability is conditional, depending on the mesh size of 
the spatial discretization and the time step of the time integration B 0]. Furthermore, in practice, the amount of 
computational work required to solve the time-dependent Maxwell equations by present FDTD techniques 0, 
||, ||, 9| , [To|, O , [l2] | prohibits applications to a class of important fields such as bioelectromagnetics and VLSI 
design jS , |13|, |14j. The basic reason for this is that the time step in the FDTD calculation has to be relatively small 
in order to maintain stability and a reasonable degree of accuracy in the time integration. Thus, the search for new 
algorithms that solve the Maxwell equation focuses on removing the conditional stability of FDTD methods and on 
improving the accuracy/efficiency of the algorithms. 

A systematic approach to construct unconditionally stable algorithms is to employ a Suzuki product-formula |l5| ] 
to approximate the time evolution operator |jq| . In the case of EM fields, the latter is the matrix exponential of a 
skew-symmetric matrix and the approximations take the form of products of orthogonal transformations [ fill |l2]| . The 
resulting numerical algorithms are unconditionally stable by construction Jl6| . |l7| |. 

The spectral-domain split-operator technique proposed in Ref. [ flO| is one of the many forms that results from the use 
of the Lie- Trotter-Suzuki product formulas JT~5[ | . This technique makes use of Fast Fourier Transforms to compute the 
matrix exponentials of the displacement operators. The choice made in Ref. jl^] yields an approximation to the time- 
evolution operator that is no longer orthogonal and hence unconditional stability is not automatically guaranteed [|l8| . 
In contrast, the methodology proposed in Refs. |n], 12 yields efficient, explicit, unconditionally stable schemes that 
operate on the EM fields defined on the real space grid only. These algorithms naturally allow for the spatial variations 
in both the permittivity and the permeability. 

The Suzuki product-formula approach also provides a unified framework to construct and analyse other time 
stepping algorithms |l6l |l9) . To illustrate this point we show that the conventional Yee algorithm and the alternating- 
direction-implicit (ADI) time-stepping algorithms || [7J |s|, [l9| fit into this framework. Furthermore we propose new 
variants of the Yee algorithm. 
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Another route to improve upon the accuracy /efficiency of time-integration schemes is to make use of the Chebyshev 



polynomial expansions of the matrix exponential |2d, 
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In this paper we also discuss a one-step 



algorithirijbased on Chebyshev polynomials, to solve the time-dependent Maxwell equations for (very) large time 
steps [f|| g7|. 

The main purpose of this paper is to review the basic ideas behind the recent developments in numerical algorithms 
to solve the time-dependent Maxwell equations and to compare the virtues and shortcomings of the different methods. 
The plan of the paper is as follows. In Section |j| we briefly discuss the basic physical symmetries of the time-dependent 
Maxwell equations. The general framework to construct time-integration algorithms is laid out in Section HI. We 
also pay attention to the numerical treatment of the current source term. In Section IV we use the simplest case of 
the time-dependent Maxwell equations to illustrate how the various algorithms can be implemented. We explicitly 
show how the conventional Yee algorithm naturally fits into this framework and, by minor modification, construct 
second-order and fourth-order time-accurate schemes that do not require the use of staggered-in-time fields, nor extra 
memory to store intermediate results. Then we recall the steps to construct the unconditionally stable algorithms 
proposed in Ref. |ll], [l^] and analyse a modification to improve the time-integration accuracy. Finally we discuss the 
implementation of the ADI and one-st ep algorithms. A discussion of the results of numerical experiments and our 
conclusions are given in Section M and VI respectively. 



II. THEORY 



We consider EM fields in linear, isotropic, nondispersive and lossless materials. The time evolution of EM fields in 
these systems is governed by the time-dependent Maxwell equations |l]] . Some important physical symmetries of the 
Maxwell equations can be made explicit by introducing the fields 

X{t) = 0*H(t) and Y(t) = \/iE(i) . (1) 

Here, H(£) = (H x (r, i), H y (r, t), H z (r, t)) T denotes the magnetic and E(<) = (E x (r, i), E y (r, t), E z (r, t)) T the electric 
field vector, while fi — fi(r) and e = e(r) denote, respectively, the permeability and the permittivity. In the absence 
of electric charges, Maxwell's curl equations Q read 

where J(t) = (J x (y, t), J y (r, i), J z (r, t)) T represents the source of the electric field and TL denotes the operator 




Writing Z(i) — (X(i), Y(t)) T it is easy to show that Tt is skew symmetric, i.e. Ti T = —H, with respect to the inner 
product (Z(f)|Z'(f)) = f v Z T (t) ■ Z'(t) dr, where V denotes the system's volume. In addition to Eq.(||), the EM fields 
also satisfy V • (^/JIX.(t)) = and V • (y/eY(t)) = 0. Throughout this paper we use dimensionless quantities: We 
measure distances in units of A and expresss time and frequency in units of A/c and c/A, respectively. 

A numerical algorithm that solves the time-dependent Maxwell equations necessarily involves some discretization 
procedure of the spatial derivatives in Eq. (g) . Ideally, this procedure should not change the basic symmetries of the 
Maxwell equations. We will not discuss the (important) technicalities of the spatial discretization (we refer the reader 
to Refs. ||, ||) as this is not essential to the discussion that follows. On a spatial grid Maxwell's curl equations (||) 
can be written in the compact form [ pd| 

^-*(t) = H*(t)-<f>(t)- (4) 

The vector vE'(i) is a representation of Z(t) on the grid. The matrix H is the discrete analogue of the operator (g), 
and the vector 4>(t) contains all the information on the current source 3(t). The formal solution of Eq. (|4|) is given by 

*(i) = U{t)*(Q) - [ U(t-u)®{u)du, (5) 
Jo 



where 
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U(t) = e tH , (6) 

denotes the time-evolution matrix. If the discretization procedure preserves the underlying symmetries of the time- 
dependent Maxwell equations then the matrix H is real and skew symmetric ||lT| , implying that U(t) is orthogonal |p8| . 
Physically, the orthogonality of U(t) implies conservation of energy p"l|. 



III. TIME INTEGRATION ALGORITHMS 



There are two, closely related, strategies to construct an algorithm for performing the time integration of the time- 
dependent Maxwell equations defined on the grid p7|. T he traditional approach is to discretize (with increasing level 
of sophistication) the derivative with respect to time Win - The other is to approximate the formally exact solution, i.e. 
the matrix exponential U (t) = e tH by some time evolution matrix U (t) jl6[ [l7] . We adopt the latter approach in this 
paper as it facilitates the construction of algorithms with specific features, such as unconditional stability |fTT[ [Tof . 

If the approximation U(t) is itself an orthogonal transformation, then ||f/(t)|| = 1 where \\X\\ denotes the 2-norm 
of a vector or matrix X j28|. In the absence of source terms (i.e. <fr(i) = 0) this implies that = ||f/(t)^(0)|| = 

||^(0)||, for an arbitrary initial condition ^(O) and for all times t and hence the time integration algorithm defined 
by U (t) is unconditionally stable by construction (l6[ [l7) . 

In the presence of current sources, for general U(t), it follows immediately from Eq.((|) that 

||*(t)|| < ||*(t)||+e^||*(0)||+^||*( W )||d^, (7) 

where ||E/(u) — f (u)|| < eTor < u < t and eTs a measure for the accuracy of the approximation U(t). 
From Eq.(||) it follows that the EM fields *f?(t) change according to 

ft+T 

¥(t + r) = e TH V{t) - J e {t+T - u)H ®{u)du. (8) 

In the time-stepping approach we approximate the source term in Eq.(|^) by the standard 3-point Gauss-Legendre 
quadrature formula [ p9[ 
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*(* + r) = e TH *(t) + \Y< me (1+x > )TH/2 $>(t + (1 + x z )t/2) + 0{t 7 ), (9) 

i=0 

where xq, x\,X2 are the zeros of the Legendre polynomial Pz{x) = x(5x 2 — 3)/2 and Wi = 8/(1 — x?)(15xf — 3) 2 [ p9| . 
In practice we replace e ( 1 + x i) TH / 2 [ n Eq.Q by an approximation U((l + Xi)r/2). 

We now consider three options to construct the approximate time evolution matrix U (t) . We exclude from the 
discussion the exceptional cases for which the matrix exponential U{t) — e tH can be calculated directly, as these 
are usually of little relevance for realistic problems. The first approach yields the conventional Yee algorithm, a 
higher-order generalization thereof, and the unconditional schemes proposed in Ref. |Tl"[ . The second option is to use 
rational approximations to the exponential, yielding the standard ADI methods. Finally, the Chebyshev polynomial 
approximation to the matrix exponential is used to construct a one-step algorithm. 

A. Suzuki product-formula approach 

A systematic approach to construct approximations to matrix exponentials is to make use of the Lie- Trotter-Suzuki 
formula |TJ 

/ v \ m 

e tH = e t(H 1+ ...+H p ) = lim fj^/m > (10) 

\i=l / 
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and generalizations thereof ||l], |3^] . Expression Eq. ([To]) suggests that 

Ui{t) = e THl ...e rH », (11) 

might be a good approximation to U(t) if r is sufficiently small. Applied to the case of interest here, if all the Hi 
are real and skew-symmetric U\(t) is orthogonal by construction and a numerical scheme based on Eq. (|ll]) will be 
unconditionally stable. For orthogonal matrices U(t) and U\{t) it can be shown that |l6| 

\\U{T)-U x (T)\\< T ^Y,\\\ H i> H Mi (12) 

i<j 

where [Hi,Hj] = HiHj — HjHi is, in general, non-zero. Relaxing the condition that U(t) and U\(t) are orthogonal 
matrices changes the r dependence in Eq. ([l2|) but for small r the error still vanishes like r 2 |52). From Eq. ( |l2| ) it 
follows that, in general, the Taylor series of U(t) and U\{t) are identical up to first order in r. We will call U\(t) a 
first-order approximation to U(t). 

The product-formula approach provides simple, systematic procedures to improve the accuracy of the approximation 
to U (r) without changing its fundamental symmetries. For example the matrix 

U 2 (t) = U 1 (-t/2) T U 1 (t/2) = e TH " /2 . . . e THl / 2 e rHl/2 . . . e rH ^ 2 , (13) 
is a second-order approximation to U(t) (3l], |32|. If U\{t) is orthogonal, so is ^(t). For orthogonal U 2 {t) we have 



\\U(t) - [U 2 (t/m)] m \\ <c 2 r 2 t, 



(14) 



where c 2 is a positive constant |iq ]. 

Suzuki's fractal decomposition approach |?2[ gives a general method to construct higher-order approximations based 
on U 2 (t) (or U\(t)). A particularly useful fourth-order approximation is given by |32j 



Ui{r) = U 2 (aT)U 2 (aT)U 2 ((l - 4a)r)C/ 2 (ar)C/ 2 (ar), 



(15) 



where a = 1/(4 — 4 1 / 3 ). The approximations Eqs.(pd|) and (Il3|), and (|T^) have proven to be very useful in many 



ipplications @ ||, || || ||, || 

solving the time-dependent Maxwell equations. As before, for orthogonal U^t) we have pi 



§ and, as we show below, turn out to be equally useful for 



\\U(t) - [[/ 4 (t/m)ni <C4tH, (16) 

where C4 is a positive constant. 

As our numerical results (see below) show, for sufficiently small r, the numerical error of a time integrator vanishes 
with t according to the r-dependence of the corresponding rigorous bound, e.g. Eqs.(^2|),(p^), or (|l6|). Our experience 
shows that if this behavior is not observed, there is a fair chance that the program contains one or more errors. 

In practice an efficient implementation of the first-order scheme is all that is needed to construct the higher-order 
algorithms Eqs.(|l^) and (|l5|). The crucial step of this approach is to choose the H.^'s such that the matrix exponentials 
exp(riii), exp(r_ff p ) can be calculated efficiently. This will turn the formal expressions for U 2 (t) and U^t) into 
efficient algorithms to solve the time-dependent Maxwell equations. 



B. ADI algorithms 



Instead of hunting for a decomposition that leads to matrix exponentials cxp(ri?i), exp(rH p ) that are easy to 
compute, one can opt for an algorithm in which each of these exponentials is calculated approximately. In principle 
this might be beneficial because there is more flexibility in decomposing H. The standard strategy, preserving 



the symmetry of Hi, H p is to use rational (Pade) approximations to the exponential |17 . For instance, the 



approximation e x sa (l + a;/2)/(l — x/2) with some decompostion H = Hi + H 2 yields the second-order- accurate ADI 
algorithm @, || g| 



U* di {t) = (I - T Hx/2)-\l + tH 2 /2)(I -tH 2 /2)-\I + tH x /2), 



(17) 
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FIG. 1: Dependence of the Bessel function J n (z = 200) on the order n. 



where / is the identity matrix. As the subscript indicates, the algorithm (17) is second-order accurate in time. For 
general skew-symmetrix H\ and H 2 , it is easy to show that the algorithm U£^^(t) is unconditionally stable. Following 
Rcf. p9| we rearrange factors and obtain 



|| [u£ DI ( T )] m \\ = \\(I-tH 1 /2)- 1 X 2 X 1 X 2 ...X 1 X 2 (I- 
< - tH 1 /2)- 1 \\\\X 2 X 1 X 2 . . . X X X 2 \ 



■tH x /2)\\ 

1(7 + ^/2)11 = -r^/2)- 1 ! 



1(1 + ^/2)1 (18) 



We used the fact that for skew-symmetric Hi, Xj = (I — tHi/2) 1 (I + rHi/2) is orthogonal and that 
\\X 2 X\X 2 ... XiX 2 \\ 2 = 1. If X is skew-symmetric, it's eigenvalues are pure imaginary and therefore (I — X)^ 1 
is non-singular. Hence, for any number of time steps m, || \U 2 DI {t)\ \\ < C, where C is some finite positive 
constant, proving that the U 2 D1 (r) is unconditionally stable in the Lax-Richtmyer sense []l7f . 

The matrix inversions appearing in Eq.(|l^) suggest that for practical purposes the implicit method U^ d ^(t) will 
not be very useful unless I — tH\ / 2 and I — tH 2 / 2 take special forms that allow efficient matrix inversion |17], [42| . 

C. One-step algorithm 

The basic idea of this approach is to make use of extremely accurate polynomial approximations to the matrix 
exponential. First we use the Chebyshev polynomial expansion to approximate U (t) and then show how to treat the 
source term in Eq. (||). We begin by "normalizing" the matrix H. The eigenvalues of the skew-symmetric matrix 
H are pure imaginary numbers. In practice H is sparse so it is easy to compute ||i7||i = maxj J+ I-Hijl- Then, by 
construction, the eigenvalues of B = —iH/\\H\\i all lie in the interval [—1,1] [§8|. Expanding the initial value *(0) 
in the (unknown) eigenvectors hj of B, we find from Eq. (||) with <&(t) = 0: 



*(t) = e lzB *(0) 
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'b,(b,l*(0)), 



(19) 



where z = £||i7||i and the bj denote the (unknown) eigenvalues of B. There is no need to know the eigenvalues 
and eigenvectors of B explicitly. We find the Chebyshev polynomial expansion of U(t) by computing the expansion 
coefficients of each of the functions e tzbj that appear in Eq. (|l9|). In particular, as — 1 < bj < 1, we can use the 
expansion [^9| e tzbj = Jq(z) + 2 J2T=i ^ Jk(z)Tk(bj) , where Jk(z) is the Bessel function of integer order k to write 
Eq. (Ill) as 



¥(t) 



J (z)I + 2j2j k (z)f k (B) 



k=l 



*(o). 



(20) 



Here Tk(B) = i k T k (B) is a matrix- valued modified Chebyshev polynomial that is defined by the recursion relations 



f Q (B)^{0) = *(0) , Ti(B)*(0) = zB*(0) 



(21) 
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and 



f fc+1 (B)*(0) - 2*BT fc (B)*(0) +f fc _ 1 (B)*(0) . 



(22) 



for fc > 1. 

As ||Tfe(B)|| < 1 by construction and | J k (z)\ < \z\ k /2 k kl for z real |^9|, the resulting error vanishes exponentially 
fast for sufficiently large K. Thus, we can obtain an accurate approximation by summing contributions in Eq. ( pc| ) 
with k < K only. The number K is fixed by requiring that | Jk{z)\ < K for all k > K . Here, K is a control parameter 
that determines the accuracy of the approximation. For fixed k, K increases linearly with z — f||-ff||i (there is no 
requirement on t being small). From numerical analysis it is known that for fixed K, the Chebyshev polynomial 
is very nearly the same polynomial as the minimax polynomial | ff2| , i.e. the polynomial of degree K that has the 
smallest maximum deviation from the true function, and is much more accurate than for instance a Taylor expansion 
of the same degree K. In Figg we show a plot of J n (z — 200) as a function of n to illustrate these points. From 
Fig.|l| it is clear that the Chebyshev polynomial expansion will only be useful if K lies to the right of the right-most 
extremum of J n (z = 200), i.e. K has to be larger than 200 in this example. 

We now turn to the treatment of the current source J(t). The contribution of the source term to the EM field at 
time t is given by the last term in Eq. (0) . For simplicity we only consider the case of a sinusoidal source 



J(r, t) = 0(i)0(T - i)s(r) sin(ta), 



(23) 



where s(r) specifies the spatial distribution and f2 = 2irf s the angular frequency of the source. The step functions 
0(i) and 0(T — t) indicate that the source is turned on at t — and is switched off at t = T. Note that Eq. ( p3|) may 
be used to compose sources with a more complicated time dependence by a Fourier sine transformation. 
The formal expression for the contribution of the sinusoidal source (B3h reads 



f e (t - u)H ${u) du = (n 2 + # 2 )- V- T ')» x (ne T ' H - n cos nr - h sin qt')s 

Jo 

= f(H,t,T',n)S, 



(24) 



where T' = min(i, T) and $(zt) = 0(t)0(T — t)sin(Ot)S. The vector S represents the spatial (time-independent) 
distribution s(r) and has the same dimension as *(0). The coefficients of the Chebyshev polynomial expansion of the 
formal solution (|24| ) are calculated as follows. First we repeat the scaling procedure described above and substitute in 
Eq. @ H = iaf||-ff||i, t = T' = Z'/\\H\\i, and Q = w\\H\\i. Then, we compute the (Fast) Fourier Transform 

with respect to x of the function f(x,z,Z',w) (which is non-singular on the interval —1 < x < 1). By construction, 
the Fourier coefficients Sk(t\\H ||i) are the coefficients of the Chebyshev polynomial expansion p9[ . 

Taking into account all contributions of the source term with k smaller than K' (determined by a procedure similar 
to the one for K), the one-step algorithm to compute the EM fields at time t reads 



*(*) 



K 



J (t\\H\\x)I + 2Y,Jk{t\\H\\ 1 )f k {B) 



k=l 



*(0) 



A" 



SoitWHh)! + 2Y,S k {t\\H\\ 1 )f k {B) 



k=l 



(25) 



Note that in this one-step approach the time dependence of the source is taken into account exactly, without actually 
sampling it. 

In a strict sense, the one-step method does not yield an orthogonal approximation. However, for practical purposes 
it can be viewed as an extremely stable time-integration algorithm because it yields an approximation to the exact 
time evolution operator U(t) = e tH that is exact to nearly machine precision, i.e. in practice the value of e m Eq.(j7j) 
is very small. This also implies that within the same precision V • (fiR(t)) = V • (fiR{t = 0)) and V • (eE(i)) = 
V • (eE(t = 0)) holds for all times, implying that the numerical scheme will not produce artificial charges during the 
time integration 0, E| . 



IV. IMPLEMENTATION 



The basic steps in the construction of the product-formula and one-step algorithms are best illustrated by considering 
the simplest case, i.e. the Maxwell equations of a ID homogeneous problem. From a conceptual point of view nothing 
is lost by doing this: the extension to 2D and 3D nonhomogeneous problems is straigthforward, albeit technically 
non-trivial 0, g g |7|. 

We consider a system, infinitely large in the y and z direction, for which e — 1 and fi = 1. Under these conditions, 
the Maxwell equations reduce to two independent sets of first-order differential equations H, the transverse electric 
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E z 



n-2 n-1 



FIG. 2: Positions of the two TM-mode field components on the one-dimensional grid. The distance between two next-nearest 
neighbors is denoted by 8. 



(TE) mode and the transverse magnetic (TM) mode [Q. As the equations of the TE- and TM-mode differ by a sign 
we can restrict our considerations to the TM-mode only. The magnetic field H y (x, t) and the electric field E z (x, 1) of 
the TM-mode in the ID cavity of length L are solutions of 



d d 
—H y {x,t) = —E„(x,t), 



i) 



d 



gt E z (x,t) = —H v (x,t)-J z (x,t), 



(26) 
(27) 



subject to the boundary condition E z (0, t) = E z (L,t) = 00. Note that the divergence of both fields are trivially 



zero. 



Following Yee J4|, to discretize Eqs.([26|) and (|27|), it is convenient to assign H y to odd and E z to even numbered 
lattice sites, as shown in Fig. |^. Using the second-order central-difference approximation to the first derivative with 
respect to a;, we obtain 



0_ 

Of 



H y (2i + l,t) = S-\E z (2i + 2,t) -E g (2i,t)), 



d_ 

df 



E z {2t,t) = 5- 1 (H y {2t + l,t)- H y (2i- l,t))- J z (2i,t), 



(28) 
(29) 



where we have introduced the notation A(i,t) = A(x = iS/2,t). The integer i labels the grid points and S denotes 
the distance between two next-nearest neighbors on the lattice (hence the absence of a factor two in the nominator) . 
We define the n-dimensional vector *f?(t) by 



*(*,*) 



__ J H y (i, t), i odd 
E z (i,t), i even 



(30) 



The vector \&(i) contains both the magnetic and the electric field on the lattice points i = 1, . . . , n. The i-th element of 
*f?(t) is given by the inner product t) = e[ ■ *S?(t) where denotes the i-th unit vector in the n-dimensional vector 
space. Using this notation (which proves most useful for the case of 2D and 3D for which it is rather cumbersome to 
write down explicit matrix representations), it is easy to show that Eqs.(28) and ( p9|) can be written in the form (^) 
where the matrix H is given by 



71—1 



/ 5- 1 

-5- 1 



e i e i+l e 



i=l 



s- 1 

-5- 1 



\ 



5- 1 
-5- 1 / 



(31) 



We immediately see that H is sparse and skew-symmetric by construction. 



A. Yee- type algorithms 

First we demonstrate that the Yee algorithm fits into the product-formula approach. For the ID model ([n]) it is 
easy to see that one time-step with the Yee algorithm corresponds to the operation 



UY ee {r) = (I + tA){I -tA t ) = e TA e~ TAT , 



(32) 
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where 



.4 



-i+l, 





s- 1 





























-5- 1 







































V : 











(33) 



and we used the arrangements of H and E fields as defined by Eq. ( ]30| ) . We use the notation JTf to indicate that the 
stride of the summation index is two. Note that since A 2 = we have e rA = 1 + tA exactly. Therefore we recover the 
time-step operator of the Yee algorithm using the first-order product formula approximation to e rH and decomposing 
H = A — A T . However, the Yee algorithm is second-order, not first order, accurate in time J2|, ||. This is due to the 
use of a staggered grid in time • To perform one time step with the Yee algorithm we need to know the values 
of E z (t) and H y (t + r/2), not H y (t). Another method has to supply the f/y-field at a time shifted by r/2. 

Within the spirit of this approach, we can easily eliminate the staggered-in-time grid at virtually no extra com- 
putational cost or progamming effort (if a conventional Yee code is available) by using the second-order product 
formula 



Ul ee {r) = e TA / 2 e- rAT e rA ' 2 = (I + tA/2){I - tA t )(I + tA/2). (34) 

The effect of the last factor is to propagate the H y -He\d by t/2. The middle factor propagates the -E^-field by r. 
The first factor again propagates the H y field by r/2. In this scheme all EM fields are to be taken at the same 
time. The algorithm defined by Urf ee (r) is second-order accurate in time by construction En]. Note that e rA ^ 2 is 
not orthogonal so nothing has been gained in terms of stability. Since (J7^ ee (r)) m = e~ rA l 2 (J7^ ee (r)) m e +Tj4 / 2 , we 
see that, compared to the original Yee algorithm, the extra computational work is proportional to (I + 2/m), hence 
negligible if the number of time steps m is large. 



According to the general theory outlined in Sec. Ill, the expression 



Ul ee (r) = U^(ar)U^(ar)U^%(l - Aa)r)Uf (ar)U^ (or) , (35) 

defines a fourth-order accurate Yee-like scheme, the realization of which requires almost no effort once UY ee has 
been implemented. It is easy to see that the above construction of the Yee-like algorithms holds for the much more 
complicated 2D, and 3D inhomogeneous case as well. Also note that the fourth-order Yee algorithm Uj ee does not 
require extra storage to hold field values at intermediate times. 



B. Unconditionally stable algorithms 



Guided by previous work on Schrodinger and diffusion problems |T^, |3^, ^] , we split H into two parts 
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such that H = Hi + H 2 - In other words we divide the lattice into odd and even numbered cells. Clearly both Hi and 
H 2 are skew-symmetric block-diagonal matrices, containing one lxl matrix and (n — l)/2 real, 2x2 skew-symmetric 
matrices. According to the general theory given above, the first-order algorithm is given by 



Ui(t) 



e rH le rH 2 



• J exp [rS 1 (e t e 




Y\_ exp [rS 1 (e 



i+l c i+2 



0] 



(38) 



To derive Eq.(^8|) we used the block-diagonal structure of Hi and 7?2 (see Eqs.(|3^) and ([37])) and obtained an exact 
expression for U\(t) in terms of an ordered product of matrix exponentials: the order of the matrix exponentials 
between each pair of curly brackets is irrelevant as these matrices commute with each other. Each of these matrix 
exponentials only operates on a pair of elements of S&(t) and leaves other elements intact. The indices of each of these 
pairs are given by the subscripts of e and e T . From Eq.(|38|) it is clear what a program should do: Make loops over 
i with stride 2. For each i pick a pair of elements from according to the subscripts of e and e T , compute (or 
fetch from memory) the elements of the plane rotation (see Eq. ([39])), perform the plane rotation, i.e. multiply the 
2x2 matrices and the vectors of length two, and overwrite the same two elements. As the matrix exponential of a 
block-diagonal matrix is equal to the block-diagonal matrix of the matrix exponentials of the individual blocks, the 
numerical calculation of e rHl (or e rH2 ) reduces to the calculation of (n — l)/2 matrix exponentials of 2 x 2 matrices. 
The matrix exponential of a typical 2x2 matrix appearing in e rHl or e rH2 is given by 



exp 




*(*,*) 



cos a 
— sin a 



sin a 

cos a 



(39) 



Using the algorithm to compute Eq.(|38|), it is easy to construct the unconditionally stable, higher-order algorithms 
U 2 {t) and U a {t), see Eq.@ and EqJ@. 

Obviously, the decomposition into Hi Eq.(|3^) and H 2 Eq.(|3^) yields the most simple real-space algorithm. It is 
not difficult to imagine that a better but slightly more complicated algorithm can be constructed by using blocks of 
3x3 instead of 2 x 2 matrices. Thus we are lead to consider the decomposition 



n-2 



#3 = S ( e » e f+l + e i+l e f+2 - ' 



i e„ - e„ 



/ 5- 1 
1 -5- 1 
-5- 1 






V 



n— 4 



■'i+3 c i+4 



/0 





-5- 1 










s- 1 














6- 1 


-5- 1 




V 









-5- 1 







s- 1 






\ 



(40) 



(41) 



where the double prime indicated that the stride of the index i is three. Obviously both H% and H A are skew-symmetric 
block-diagonal matrices, build from the 3x3 skew-symmetric matrix 



(42) 




As B 3 = -2B we have 



1 + sB + cB 2 




(43) 



where s = sin(v / 2r) and c = sin 2 (r/\/2). In practice, using the 3x3 instead of the 2x2 decomposition is marginally 
more difficult. We will denote the corresponding second and fourth-order algorithm by t/| x3 and f/f x3 respectively. 
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FIG. 3: The field H y (x, t = 100) generated by a current source at x = 125 that oscillates at frequency / s = 1 during the interval 
< t < 6, as obtained by the one-step algorithm with K = 2103 [K = in this case). 



C. ADI algorithm 

For the tri-diagonal matrix (|3l|), the ADI algorithm reduces to the Cranck-Nicholson method p2|| . The tri-diagonal 
structure of the matrix H permits the calculation of (I — TH/2)^ li Sf in 0(n) operations by standard linear algebra 
methods f!§. 

D. One-step algorithm 

The one-step algorithm is based on the recursion (see Eqs.(]2l]) and (p2[)) 

2H 

= TEnr^k + *fc-i- (44) 

ll-nlli 

Thus, the explicit form Eq.(pT|) is all we need to implement the matrix- vector operation (i.e. <— H^) that enters 
Eq.@. 

The coefficients Jh{z) and Sk(t) (see Eq.(p5|)) should be calculated to high precision. Using the recursion relation 
of the Bessel functions, all K coefficients can be obtained with O(K) arithmetic operations p2| |. The numbers Sk(t) 
can be calculated in 0(K log K) by standard Fast Fourier transformation techniques. Clearly both computations are 
a neglible fraction of the total computational cost for solving the Maxwell equations. 

Performing one time step amounts to repeatedly using recursion ( p2f ) to obtain Tk(B)^(0) for k = 2, . . . , K, multiply 
the elements of this vector by Jh{z) (or Sk{z)) and add all contributions. This procedure requires storage for two 
vectors of the same length as ^(0) and some code to multiply such a vector by the sparse matrix H. The result 
of performing one time step yields the solution at time t, hence the name one-step algorithm. In contrast to what 
Eqs. ( pl| ) and (^2|) might suggest, the algorithm does not require the use of complex arithmetic. 

In the sequel, the caret on top of a symbol indicates that the results have been obtained by means of the one-step 
algorithm. 

V. NUMERICAL EXPERIMENTS 

Except for the conventional Yee algorithm, all algorithms discussed in this paper operate on the vector of fields 
defined at the same time t. We use the one-step algorithm (with a time step r/2) to compute E z (t/2) and H y (r/2). 
Then we use E z (0) and H v (t/2) as the initial values for the Yee algorithm. In the presence of a current source, there 
are some ambiguities with this procedure as it is not obvious how to treat the source term in Eq.(^). In order to 
permit a comparison of the final result of the conventional Yee algorithm with those of the other methods, we use 
the one-step algorithm once more to shift the time of the H y field by —r/2. This procedure to prepare the initial 
and to analyse the final state of the Yee algorithm does in fact make the results of the Yee algorithm look a little 
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TABLE I: The error ||*(t) - *(t)||/||*(t)|| at time t = 100 as a function of the time step r for eight different FDTD algorithms. 
The current source is positioned at the center of the system (see Figj^), and oscillates at frequency f s = 1 during the interval 
< t < 6, *fr(t) is the vector obtained by the one-step algorithm with re = 10 -9 , using K' — 2103 matrix-vector operations 
\E"' <— M^Sf. Yee: \&(t) obtained by the Yee algorithm g, Other columns: *ff(t) obtained by the algorithms indicated. 
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FIG. 4: The data presented in Table |] plotted on a double logarithmic scale. Lines are guide to the eye. 



more accurate than they would be if the exact data of the r/2-shifted fields are not available. Thus, the results on 
the errors of the conventional Yee algorithm presented in this paper give a too optimistic view on the accuracy of 
this algorithm but we nevertheless adopt the above procedure to make a quantitative comparison between the various 
algorithms. 

We define the error of the solution ^(t) for the wave form by ||^(t) — *&(i)||/||\&(f)|| where *&(t) is the vector of 
EM fields obtained by the one-step algorithm. Thereby we have already assumed that the one-step algorithm yields 
the exact (within numerical precision) results but this has to be demonstrated of course. A comparison of the results 
of an unconditionally stable algorithm, e.g. U4 with those of the one-step algorithm is sufficient to show that within 
rounding errors the latter yields the exact answer. Using the triangle inequality 



||*(t)-*(t)||<||*(i)-*(*)|| + ||*(i)-*(*)||, 



(45) 
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TABLE II: The error ||*(i)-*(t)|l/||*(t)|| at time t = 100 as a function of the time step r for eight different FDTD algorithms. 
The system is the same as in Fig^| and Table 0. The initial values of the EM fields are random, distributed uniformly over 
the interval [-1,1]. *&{t) is the vector obtained by the one-step algorithm k, — 10~ 9 , using K = 2080 matrix-vector operations 
«— M^Sf. Yee: ^f(t) obtained by the Yee algorithm j^, ^, Other columns: 4*(t) obtained by the algorithms indicated. 
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FIG. 5: The data presented in Table [o] plotted on a double logarithmic scale. Lines are guide to the eye. 



and the rigorous bound 



\\*{t) - *(*)|| < c 4 tH (j|*(0)|| + J* ||J(u)||du) , (46) 

we can be confident that the one-step algorithm yields the numerically exact answer if i) Eq.djpj) is not violated and 
ii) if H* ft) - *(i)|| vanishes like t 4 . 

In Fig.0 we show a typical result of a one-step calculation on a grid of n = 5001 sites with S = 0.1 (corresponding to 
a physical length of 250.05), and a current source placed at i = 2500 to eliminate possible artifacts of the boundaries. 
The frequency of the source is set to one (/ s = 1) and the number of periods the source radiates is set to six (i.e. 
T = 6). In Table | (Fig.||) we present results for the errors, as obtained by repeating the simulation shown in Fig.|| 
using eight different FDTD methods. In Tables |l| and III (Figs.|| and || respectively) we shown similar results but 
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TABLE III: The error ||*(t) - *(i)||/||*(i)|| at time t = 100 as a function of the time step r for eight different FDTD 
algorithms The system is the same as in Fig^| and Table [l| The initial state of the EM fields is a Gaussian wave packet 
(E z (t) = exp(— (x — xq — t) 2 /a 2 ) with a width a = 4, and its center xo = 125 positioned at the middle of the system (see Pig.^j) . 
9(t) is the vector obtained by the one-step algorithm with k = 10~ 9 , using K = 2080 matrix-vector operations 9' <— M3". 
Yee: 9(t) obtained by the Yee algorithm 0,0,0]; Other columns: 9(t) obtained by the algorithms indicated. 
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instead of using a current source, a random wave form (Table || and Gaussian wave packet (Table |ffl|) 

was taken as 

the initial condition. 

From the data in Tables |, || and III we conclude that the error of algorithm U4 vanishes like r 4 , demonstrating that 
the one-step algorithm yields the numerically exact result (see Eqs.[l5]and [l6|). The results presented in Tables || and 
[II have been obtained by using a vector of initial values that is normalized to one, i.e. ||\lf(0)|| = 1. As ||ifr(t)|| = 1 
for < t < 100 to at least 10 digits, - *(t)||/||*(f)|| = - 4f(t)\\ for all entries in Tables || and [fij. The 

high precision of the one-step algorithm also allows us to use it for genuine time stepping with arbitrarily large time 
steps, this in spite of the fact that strictly speaking, the one-step algorithm is not unconditionally stable. 

The data in Tables |, |l| and III suggests that there does not seem to be a significant difference between the 
conventional Yee algorithm and its variant UY ee but in fact there is. The time evolution matrix corresponding to 
the Yee and the l/Y ee algorithm is not orthogonal. Therefore the energy of the electromagnetic field (|| ^(t)!! 2 ) is not 
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FIG. 7: The energy of the EM fields W = * T (t) ■ ¥(t) as the function of time as obtained by the Yee (solid line), U¥ ee 
(dashed line), and U4 (dotted line) algorithm for a ID cavity of size 48.05 (n = 97 mesh points), a mesh size 5 — 0.1 and a time 
step r = 0.01. 



conserved. Furthermore, in the conventional Yee algorithm, the E and //-fields are time-shifted by r/2 with respect 
to each other. These artifacts of are less prominent if we use lf% ee instead of the Yee algorithm. In Fig.f?] we show 
results of the time evolution of the total energy of the EM field, for a system of n=97 sites, a mesh size 8 = 0.1 and 
a time step of t = 0.01. A normalized (||\&(0)|| 2 = 1) random initial condition was used. Furthermore, for this type 
of application, it makes no sense to invoke the procedure described at the beginning of this section to time-shift one 
of the EM-fields by r/2: as this operation has to be performed at each time step and is computationally expensive 
(because it is numerically exact), we could as well use the same numerically exact procedure for time stepping itself. 

For the Yee algorithm, the fluctuations of the energy are a factor often larger than in the case of the U^ ee algorithm. 
As expected on theoretical grounds, the U4 algorithm (dotted, horizontal line) exactly conserves the energy. The fact 
that U£ ee conserves EM-field energy much better than the Yee algorithm also has a considerable impact on the quality 
of the eigenmode distribution. The latter is obtained by Fourier transformation of ^f T (t) ■ ^(0) (see ref.|ll|] for more 
details). In Fig.|| we show the low- frequency part of the eigenmode distribution of the same system as the one of 
Fig.^. It is obvious that there is a significant improvement in the quality of the spectrum if we use C/^ ee instead of 
the Yee algorithm but for this application U4 performs much better than the Yee-type algorithms. 

Table | suggests that U2 is the least efficient of the five FDTD methods: It uses more arithmetic operations than 
the Yee algorithm and yields errors that are larger than those of the Yee algorithm. However, this conclusion is biased 
by the choice of the model problem and does not generalize. If the initial EM field distribution is random then, for 
sufficiently small r, algorithm U2 is more accurate than the two second-order accurate Yee algorithms, as is clear from 
the data in Table [0] p3fl . Also in this case, for the largest r in Table ||, the Yee and [/^ ee algorithm are operating 
at the point of instability, signaled by the fact that the norm of S&(t) grows rapidly, resulting in errors that are very 
large. From Tables [| and || one might conclude that the decomposition that generates Yee-type algorithms yields 
the least accurate approximations to the time evolution operator, although the difference is not really significant, 
but, as Table III shows, this conclusion would be wrong. If the initial state is a Gaussian wave packet that is fairly 



broad, the Yee-type algorithms arc muc h m ore accurate than the unconditionally stable algorithms employed in this 
paper. From the data in Tables |, || and III we conclude that there is no good reason to use the ADI algorithm (even 



disregarding the fact that it is slower than the other second-order methods). In general the t/| x3 (Z7f x3 ) algorithm 
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FIG. 8: The eigenvalues distribution of the matrix H, as obtained by Fourier transformation of *f? T {t) ■ ^(0), for the same 
system as in Figj^. The function VE fT (i) • *(0) is sampled at time intervals of 0.1, the total number of samples being 4096. Solid 
line: Yee algorithm; dashed line: U~2 ee algorithm; dotted line: energy conserving algorithm U4. 



performs a little better than U2 (C/4) but the gai n is marginal. In contrast to the numerical data presented in Ref.|19|, 
for all algorithms the data of Tables |,pl| and |l| clearly agree with the theoretically expected behavior of the error as 
a function of r if r is small enough K4L 

Usually if a current source is present we have >P(0) = 0. Then the one-step algorithm requires K' (sparse) matrix- 
vector operations ( 1 4 ,/ <— M\&) to compute ^(t). For a ID system the standard Yee, ~U^ ee and Uj ee ,1/2 , and U4 
algorithms perform (in worst case, without additional optimization), respectively, 1, 3/2, 6, 3/2, and 6 M \&-operations 
per time step. The one-step algorithm carries out K' — 2103 matrix-vector operations <— to complete this 
simulation. This implies that for all r < t/K', the FDTD algorithms will perform more \J r ' <— M\P operations than 
the one-step algorithm. For the data presented in this paper, this is the case if r < 0.05 for the Yee algorithm and is 
always the case for U4 because the latter uses a factor of 6 more <— M\& operations than the Yee algorithm. 



VI. CONCLUSION 



The answer to the question which of the algorithms is the most efficient crucially one depends on the accuracy that 
one finds acceptable. Taking the data of Table | as an example we see that if one is satisfied with an error of more 
than 2%, one could use the Yee algorithm. With r = 0.05 it needs 2000 time steps to find the solution at t = 100, 
close to the K' = 2103. Nevertheless we recommend to use the one-step algorithm because then the time-integration 
error is neglegible. The Yee algorithm is no competition for U4 if one requires an error of less than 1% but then U4 is 
not nearly as efficient (by a factor of about 6) as the one-step algorithm. Increasing the dimensionality of the problem 
favors the one-step algorithm [p6| , p7| . These conclusions seem to be quite general and are in concert with numerical 
experiments on ID, 2D and 3D systems ^7|. A simple theoretical analysis of the r dependence of the error shows 
that the one-step algorithm is more efficient than any other FDTD method if we are interested in the EM fields at a 
particular (large) time only p6| , p7| . This may open possibilities to solve problems in computational electrodynamics 
that are currently intractable. The Yee-like algorithms do not conserve the energy of the EM fields and therefore 
they are less suited for the calculation of the eigenvalue distributions (density of states), a problem for which the U4 
algorithm may be the most efficient of all the algorithms covered in the paper. 
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The main limitation of the one-step algorithm lies in its mathematical justification. The Chebyshev approach 
requires that H is diagonalizablc and that its eigenvalues are real or pure imaginary. The effect of relaxing these 
conditions on the applicability of the Chebyshev approach is left for future research. 

In this paper we have focused entirely on the accuracy of the time integration algorithms, using the most simple 
discretization of the spatial derivatives. In practice it is straightforward, though technically non-trivial, to treat more 
sophisticated discretization schemes || [l^] by the methodology reviewed is this paper. 



Acknowledgments 

H.D.R. and K.M. are grateful to T. Iitaka for drawing our attention to the potential of the Chebyshev method and 
for illuminating discussions. This work is partially supported by the Dutch 'Stichting Nationale Computer Faciliteiten' 
(NCF), and the EC 1ST project CANVAD. 



M. Born and E. Wolf, Principles of Optics, (Pergamon, Oxford, 1964). 

A. Taflove and S.C. Hagness, Computational Electrodynamics - The Finite-Difference Time-Domain Method, (Artech 
House, Boston, 2000). 

K.S. Kunz and R.J. Luebbers, Finite- Difference Time-Domain Method for Electromagnetics, (CRC Press, 1993). 
K.S . Yee, IEEE Transact ions on Antennas and Propagation 14, 302 (1966). 
See |ittp:/ /www. fdtd.org 



F. Zheng, Z. Chen, and J. Zhang, IEEE Trans. Microwave Theory and Techniques 48, 1550 (2000). 
T. Namiki, IEEE Trans. Microwave Theory and Techniques 48, 1743 (2001). 
F. Zheng and Z. Chen, IEEE Trans. Microwave Theory and Techniques 49, 1006 (2001). 
S.G. Garcia, T.W. Lee, S.C. Hagness, IEEE Trans, and Wireless Prop. Lett. 1, 31 (2002). 
W. Harshawardhan, Q. Su, and R. Grobe, Phys. Rev. E 62, 8705 (2000). 
J.S. Kole, M.T. Figge and H. De Raedt, Phys. Rev. E 64, 066705 (2001). 
J.S. Kole, M.T. Figge and H. De Raedt, Phys. Rev. E 65, 066705 (2002). 

O.P. Gandi, Advances in Computational Electrodynamics - The Finite- Difference Time-Domain Method, A. Taflove, Ed., 
(Artech House, Boston, 1998). 

B. Houshmand, T. Itoh, and M. Piket-May, Advances in Computational Electrodynamics - The Finite-Difference Time- 
Domain Method, A. Taflove, Ed., (Artech House, Boston, 1998). 
M. Suzuki, S. Miyashita, and A. Kuroda, Prog. Theor. Phys. 58, 1377 (1977). 
H. De Raedt, Comp. Phys. Rep. 7, 1 (1987). 

GD. Smith, Numerical solution of partial differential equations, (Clarendon Press, Oxford, 1985). 

This is due to the specific split-up adopted in [B~o|l and not an intrinsic property of the spectral-domain approach. 



R. Horvath, preprint 2002; http://www.win.tue.nl/analysis/preprints/2002.htm] 



H. Tal-Ezer, SIAM J. Numer. Anal. 23, 11 (1986). 
H. Tal-Ezer and R. Kosloff, J. Chem. Phys. 81, 3967 (1984). 

C. Leforestier, R.H. Bisseling, C. Cerjan, M.D. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, 

H.-D. Meyer, N. Lipkin, O. Roncero, and R. Kosloff, J. Comp. Phys. 94, 59 (1991). 

T. Iitaka, S. Nomura, H. Hirayama, X. Zhao, Y. Aoyagi, and T. Sugano, Phys. Rev. E 56, 1222 (1997). 

R.N. Silver and H. Roder, Phys. Rev. E 56, 4822 (1997). 

Y.L. Loh, S.N. Taraskin, and SR. Elliot, Phys. Rev. Lett. 84, 2290 (2000); ibid. Phys.Rev.Lett. 84, 5028 (2000). 

H. De Raedt, K. Michielsen, J.S. Kole, and M.T. Figge, in Computer Simulation Studies in Condensed-Matter Physics 

XV, eds. D.P. Landau et al., Springer Proceedings in Physics, (Springer, Berlin in press). 

H. De Raedt, K. Michielsen, J.S. Kole, and M.T. Figge, submitted to IEEE Antennas and Propagation, LANL prepint: 



physics/0208060 



J.H. Wilkinson, The Algebraic Eigenvalue Problem, (Clarendon Press, Oxford, 1965). 

M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, (Dover, New York, 1964). 

H.F. Trotter, Proc. Am. Math. Soc. 10, 545 (1959). 

H. De Raedt and B. De Raedt, Phys. Rev. A 28, 3575 (1983). 

M. Suzuki, J. Math. Phys. 26, 601 (1985); ibid 32 400 (1991). 

A.J. Chorin, T.J.R. Hughes, M.F. McCracken, and J.E. Marsden, Comm. Pure Appl. Math. XXXI, 205 (1978). 

H. Kobayashi, N. Hatano, and M. Suzuki, Physica A 211, 234 (1994). 

H. De Raedt, K. Michielsen, Comp. in Phys. 8, 600 (1994). 

A. Rouhi, J. Wright, Computers in Physics 9, 554 (1995). 

BA. Shadwick and W.F. Buell, Phys. Rev. Lett. 79, 5189 (1997). 

M. Krech, A. Bunker, and D.P. Landau, Comp. Phys. Comm. Ill, 1 (1998). 

P. Tran, Phys. Rev. E 58, 8049 (1998). 



17 



[40] K. Michielsen, H. De Raedt, J. Przeslawski, and N. Garcia, Phys. Rep. 304, 89 (1998). 

[41] H. De Raedt, A.H. Hams, K. Michielsen, and K. De Raedt, Comp. Phys. Comm. 132, 1 (2000). 

[42] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, Numerical Recipes, (Cambridge, New York, 1986). 
[43] This also explains why the unconditionally stable algorithms U2 and U4 yield more accurate eigenvalue distributions than 
the Yee algorithm |ll| . 

[44] The data for small r, obtained from the Yee algorithm for the case of a current source is an exception. This is due to the 
ambiguities mentioned in Sec.M. 



